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Abstract 

As an approach to describe the long-range properties of non-Abelian gauge theories at non-zero 
temperature T < T c , we consider a non-interacting ensemble of dyons (magnetic monopoles) 
with non-trivial holonomy. We show analytically, that the quark-antiquark free energy from 
the Polyakov loop correlator grows linearly with the distance, and how the string tension scales 
with the dyon density. In numerical treatments, the long-range tails of the dyon fields cause 
severe finite-volume effects. Therefore, we demonstrate the application of Ewald's summation 
method to this system. Finite-volume effects are shown to be under control, which is a crucial 
requirement for numerical studies of interacting dyon ensembles. 



1 Introduction 



Insight into the mechanisms of the QCD vacuum is not only provided by simulations of lattice 
gauge theory - an ab initio method, whose numerical results, however, are hard to interpret - 
but also by analytical non-perturbative approaches like the semiclassical one [HE]- The latter 
relies for instance on instantons, selfdual and anti-selfdual solutions of the Euclidean Yang-Mills 
equations [3]. Instantons in M 4 are localized in space and time, but also naturally contain 
long-range fields (since the Yang-Mills Lagrangian is scale-invariant): the gauge potential 
decays like the inverse of the four-dimensional distance to their center or its third power, in the 
regular and singular gauges, respectively. Semiclassically motivated models of the QCD vacuum 
based on instantons are suitable to describe certain non-perturbative effects like chiral symmetry 
breaking, but so far cannot explain confinement. For more details of instanton models we refer 
to the reviews [U [5] . 

When studying instantons or similar long-range (or "infrared") objects in a finite- volume ap- 
proximation - an unavoidable restriction for virtually every numerical approach - one expects 
severe effects: interactions with objects outside the finite volume (and their contribution to ob- 
servables) are neglected, which can introduce considerable systematic deviations from analogous 
systems with infinite extent. 

The purpose of our work is two-fold. On the one hand, we investigate confinement in a semi- 
classical approach at non-zero temperature. Guided by the invention of KvBLL-calorons with 
non-trivial holonomy [6j O E] our basic objects are dyons - the constituents of calorons. We as- 
sume maximally non-trivial holonomy in order to describe the confinement phase of the model. 
Dyons will be analytically shown to provide a confining Polyakov loop correlator already within 
the simplest non-interacting model for the low-temperature phase. 

Concerning the long-range nature, dyons are as difficult to simulate as instantons. Therefore, 
as the second part, we provide the proof-of-concept for a method capable to control the finite- 
volume effects in such systems in an efficient way: Ewald's summation method [9]. This method 
was originally developed for Coulomb interactions typical, for example, in plasma or soft matter 
physics. When using Ewald's method, the infinite space is mimicked by infinitely many replicas 
of one so-called "supercell" that contains a (for numerical simulations) feasible number of ob- 
jects/charges. Typical observables like potentials are sums, which can be split into a short-range 
and a long-range part. After rewriting the long-range sum by means of a Fourier transform, 
both sums can efficiently be computed (see Section [4.11 for details). In order to come back to 
the original system in the infinite space, only the volume of the supercell has to be extrapolated 
to infinity at the end of a computation. In Ewald's method this is a well-controlled limit in 
contrast to simpler approaches. 

In this work we apply this method to the simplest dyon model, for which a comparison with 
analytical results can be made. The advantages of the Ewald method will be essential for later 
numerical studies of interacting dyon ensembles. 

The paper is organized as follows. In Section[2]the crucial features of dyons are introduced. The 
Polyakov loop correlator in a non-interacting dyon model is analytically evaluated in Section [3l 
both in infinite and finite volume. Contact to lattice simulations is made and consequences 
for dyon models are discussed. Section [J] introduces Ewald's summation method. In Section [5] 
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numerical results are presented. Section [6] summarizes this work and opens a view to simulations 
of interacting dyon models. In Appendix [A] some integrals required in the analytical approach are 
computed, whereas Appendix [B] compares Ewald's method of summing over an infinite number 
of copies of the supercell with the result of the converging series of sums over a finite number of 
cells. 

2 Dyon gas model for SU(2) Yang-Mills theory 

The notion of finite temperature instanton solutions, traditionally called calorons jlOj . has been 
radically extended when new caloron solutions were found by Kraan and van Baal [6j [7] as well as 
Lee and Lu |5]. They consist of magnetic monopoles as constituents. The latter also carry (the 
Euclidean analog of) electric charge and will therefore be called dyons. The asymptotic Polyakov 
loop of these solutions, the trace of the so-called holonomy, is an additional external parameter 
that governs for instance how the instanton (caloron) action is shared by the constituent dyons. 

Dyons as selfdual objects at finite temperature can be obtained by considering the gauge field of 
a caloron in the limit of infinite dyon separation [6] . The dyon constituents can be understood 
as BPS monopoles interpreting the scalar Higgs field as a temporal gauge field. In the far-field 
limit, when the distance to the dyon center is large, the gauge field is Abelian along the direction 
of the asymptotic Polyakov loop ( "the color direction of the Higgs field" ) , which we take diagonal 



with T the temperature and 03 = diag(+l, —1) the third Pauli matrix. The parameter u 
specifies the holonomy. Maximally non-trivial holonomy refers to uj = 1/4 and -P(r) — > and 
is conjectured to be valid in the confined phase, where (P) = (as a quantum spatial average), 
in contrast to trivial holonomy P(r) — > ±1 valid deep in the deconfined phase. The viability 
of confinement has been shown semi-analytically, even without complete decomposition into 
constituents |llj . Further investigations supporting the conjecture have been focused on the 
quantum amplitude |12| , moduli space metric [13|, [T4"] and the vortex content of calorons |15| . 

Dyons are genuine non-Abelian objects, whose field components color-perpendicular to the 
asymptotic Polyakov loop decay exponentially (like e.g. fields of massive bosons color perpen- 
dicular to the Higgs vacuum expectation value) outside a region of size j3 = 1/T. The dyons' 
long-range fields are Abelian in the same color direction and Coulomb-like (in addition to the 
constant of Eq. ([!])): 



where r = (x,y,z) and r = |r| is the three-dimensional distance to the dyon center. With the 
help of 't Hooft's symbol one can write in a compact way 



A -> 2-kojT 03 






(3) 




(4) 
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The vector potential in this limit results in electric and magnetic fields 



e 



qr 



b = 



q (4 + 4ir6(x)5(y)Q(z)e z 



(5) 



The possible charges are q = +1 for dyons and q = — 1 for anti-dyons. The Dirac string 
singularities along the positive z- axis are artefacts of the Abelian limit. They do not need to 
concern us here. 

So far we have considered selfdual dyons, whose electric and magnetic charges are coupled as 
e = b (neglecting Dirac strings). The actual semiclassical field content dominating the partition 
function should be built from selfdual and antiselfdual dyons and antidyons. For antiselfdual 
dyons and antidyons, for which e = — b, the 't Hooft symbols ff^ are replaced by r]^ u . As we will 
argue below, all that matters for our work is the Coulomb-like decay of ao away from positive 
and negative electric charges ±q, which are placed at random positions. This will apply also 
to a mixed model including antiselfdual dyons and antidyons. In other words, for the aspects 
under study selfduality and antiselfduality and the magnetic charges are irrelevant properties 
and as a consequence, our formulation respects CP-invariance. In due course, total numbers 
and densities will refer to dyons of both magnetic charges. 

The superposition of the gauge fields of 2K dyons in the Abelian limit reads 



where rf 1 and q m = — (— l) m are the positions and charges of the i-th dyon (m = 1) and antidyon 
(to = 2), respectively. 

Like the vector potentials, interactions of monopoles or dyons behave Coulomb-like |16 [ \17 \ ITS]. 
Relying on these long-range fields, Diakonov and Petrov have presented a formal solution for the 
statistical mechanics of purely (anti)selfdual dyons [19] . later extended to both selfdualities |20j . 

The assumed moduli space metric of the dyon configurations allowed for a particular analytic 
treatment in the spirit of Polyakov's monopole confinement mechanism |21j. In an attempt 
to implement a simulation for dyon gases with this interaction, however, we have noticed that 
the metric severely suffers from non-positivity [22] . which casts doubts on the validity of the 
analytical results obtained in [19] in the context of Yang-Mills theory. 

In this paper we consider dyon ensembles without moduli space metric or other interactions, 
i.e. we perform a uniform sampling of dyon positions. We will focus on maximally non-trivial 
holonomy, u = 1/4, where both dyons and antidyons possess the same topological charge of 
2u = 1 — 2u = 1/2 of an instanton unit and hence the same action, such that they do not differ 
in their classical and quantum weight. Therefore, it is natural to use the same number of dyons 
and antidyons, i.e. an electrically and magnetically neutral ensemble and denote by = 2K 
the total number of dyons and antidyons. For other values of the holonomy, say for those close 
to maximally non-trivial, the assumption of equally frequent dyons is only a first approximation, 
arguments suggesting the contrary are discussed in references [231 [23] . 

The basic parameters of our model are the holonomy uj, the 3-dimensional density of dyons p 
and the temperature T. The scale can be set by identifying the string tension a extracted from 




(6) 
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the free energy of a static quark-antiquark pair with the corresponding lattice result as explained 
in Section 13.31 

Our primary observable is the local Polyakov loop P(r) at position r (cf. Eq. ([2]))- In the 
Abelian limit the fields are static and we need to sum the holonomy and the ao-component of 
the individual dyons as follows, 



w=l/4 



( $ ( r ) 



(7) 



with the following sum over Coulomb terms 



K 2 



*w -EE 



i=l m=l 



r — r;- 



K 

£ 

i=l L 



1 i i 9| 

r — r}\ r — rf 



(8) 



As well-known, the correlator of Polyakov loops yields the free energy of a static quark-antiquark 
pair: 



Tln(p(r)P\r') 



d 



r — r 



(9) 



From the point of view of a Coulomb gas, correlators of trigonometric functions are slightly 
exotic, but for the dyon model of QCD this is the essential correlation function probing confine- 
ment. 

In simulations using a finite number of dyons and anti-dyons the positions vf 1 are restricted to a 
finite dyon sampling volume. Then contributions from dyons outside this volume to the sum in 
Eq. ([8j) are ignored. How one can control such finite-volume effects systematically, is the main 
subject of the second part of this paper. We will resort to Ewald's summation method and 
compare it to the analytic result for Polyakov loop correlators, which are presented in the next 
section. 



3 The Polyakov loop correlator in a non-interacting dyon gas 
model 

In this section we treat the non-interacting dyon ensemble analytically. In particular we show 
the Polyakov loop correlator ([9]) from random dyons to be confining and investigate finite- volume 
effects. Interacting dyon ensembles can be reformulated as scalar theories |25|, 119] , but here - 
due to the absence of interactions - the model can be solved. In this simple system we therefore 
obtain analytic formulae for the string tension, which later will be used to set the scale and as 
a benchmark for numerical methods. 
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3.1 The correlator 



Expectation values of observables O in the ensemble with K dyons of charge +1 at positions rj 
and K dyons of charge —1 at positions r? are given by: 



O 



[ n*Wo({rJ,r?})/ / n*W= / Udr}drl0({r},rn)/V 2K 

J i=l J i=l i=l 



(10) 



where V is the spatial volume in which the 2K dyons are randomly distributed. Their density 
is 



2K 

p = T 

accordingly. 



(11) 



The Polyakov loop correlator is given by a product of cosines, see (J7J, (jHj) and Q, and can be 
rewritten as 



P(r)P(r') 
where the Coulomb sums 



1 / (a ®+ 

— ( COS iTTUJ H 

2 \ V 2T 



1 / 

H — ( cos — - — 

2 \ V2T 



$± = $(r)±$(r') = ^ 



i=l 



1 ± 1 



1 ± 1 



If 2 1 I j./ p2i 



(12) 



(13) 



contain all dyons and depend on the two measurement points r and r'. 
Rewriting 



P(r)P(r') 



1 



A-KtijJ 



2T 



1 



-e < exp I i— - ) + c.c. + -< exp 1 1— - ] ) + c.c 



2T 



(14) 



the ingredients are the following expectation values 

K 



2T 



1 



1 

V 



i=l 

ds exp 



2T V I r — r j I |r' — rj 



1 1 

± 



2TV|r-s| |r'-s| 



X c.c. 



X c.c. 



I±\ 

V 



2K 



■ (15) 



They are real and have factorized into integrals given in terms of one dyon location only: 



ds exp 



1 1 

± 



2TV|r-s| |r'-s 
The result for the Polyakov loop correlator is then 
1 



(16) 



P(r)P(r') 



cos(47ru;) 



V 



+ 



2 V V 



(17) 
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Keeping the density fixed, we can replace the number of dyons 2K and obtain 



(p{y)P(t')^ = - cos(4^u;) 
in particular at maximally non-trivial holonomy uj = 1/4 

<™>L ~ 4[(-^)T+#+^)T 

This (still exact) form with the explicit volume dependencq^ has been chosen in anticipation of 
the properties of I± discussed below. 



1 + 



V 



V 



1 

+ 2 



1 + 



V 



V 



v 



(18) 



3.2 String tension in the infinite-volume limit 

The task here will be to calculate the asymptotic behavior of the integrals I± of Eq. (|16p in 
the limit of large quark-antiquark separations. The behavior at finite separations as well as 
finite-volume corrections are investigated in the next subsection. 



By shifting and rotating the integration variable in Eq. (|16p one can see that I± are functions 
of the distance |r — r'| = d as expected. 

The integrands of both integrals /+ and J_ asymptotically approach unity, the corresponding 
(divergent) term will be canceled by V in Eq. (|19p . However, there is an important difference: 
with the relative plus sign in I + the next term in the asymptotic expansion is the monopole 
term (proportional to 2/s), while the integrand of i_ will only start with a dipole term due to 
the relative minus sign. We will show that as a consequence the first term in Eq. f)19|) vanishes 
in the infinite- volume limit, whereas the second term survives and induces the string tension. 

We consider regularized integrals in a 3-ball of radius R and fix the Polyakov loop arguments at 
r = (0, 0, +d/2) and r' = (0, 0, —d/2). Notice first that the integration variable s can be rescaled 
by the temperature 



RT 



=/± 

such that these integrals are functions of the finite-volume radius R and the separation d only, 
both in units of 1/T. In other words 

I± = ^f±(dT,RT). (21) 

In spherical coordinates the dT dependence of the distances \rT — s| and \r'T — s| becomes 
explicit: 



D±(s,9,dT) = |(0,0,±d/2)r-s| = yjs 2 T sdT cos 9 + (cff) 2 /4 . (22) 



The use of dimensionful quantities in the exponent can be avoided by normalizing with some standard volume. 
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In order to evaluate the leading terms in the radius R, we consider the i?T-derivatives of f± for 
large RT given by angle integrals on the 2-sphere s = RT: 

f± = 2ir(RT) 2 [ d6 sin# (23) 
Jo 

ex P I 7T ( Z^T~E == 5 ,? tt^ ± 



d(RT) 



2\RT v / l + d/R-cos6 + d 2 /4R 2 RT^l - d/R ■ cos 9 + d 2 /AR 2 
This can be expanded in 1/RT and to give 

^/ + = MRTf + iMRT)-2. + 0^,^) (24) 



d(RT) K ' \R 2, (RT) 2 R± 

where the first terms on the right hand sides will, of course, be the volume contributions. 
Moreover, the aforementioned difference in the two integrals concerning subleading terms is 
clearly visible. Integrating back with respect to RT then yields 

Air 

f+ = — (RT) 3 + i2ir(RT) 2 -2ir (RT) + (In RT)+g + (dT) (26) 
3 

/_ = ^(RTf + g_(dT) (27) 

where g± are i?T-independent and where we have neglected all terms vanishing as RT — > oo. 
For |i+| we finally get the following leading terms: 



4vr , R R 2 
— R 3 - 2it—r + i- 2vr— 
3 T 2 T 



V-cT 2 R, = -d T - 2 (28) 



with c and d being positive constants. In the infinite-volume limit at fixed temperature the 
contribution to the Polyakov loop correlator vanishes 

/ |/+i-y\ y / dT 2 \^ RS 

lim I + ^t = lim 1-— =- =0. (29) 

y^oo \ V J R^oo \ R 2 J 

In |i_|, on the other hand, only .R- independent terms enter the Polyakov loop correlator as 
9 



V + T 3 



V + ^ (30) 



|J_| 
and 



7 



Hence it remains to compute g_ as a function of the Polyakov loop separation d = |r — r'| (in 
units of T), which according to the above is 

g_=[ ds/expflf— L_ __ L-'lVlV (32) 



,2V|rT-s| \r'T 

The imaginary part vanishes by invar iance under reflections s — > — s. We split 



1 (2) , fres) (2) f . ( 



9 - = -^W~', ^.j^^-^-^f (33) 

such that all terms are integrable around rT and r'T. The important observation is now that 
the second-order contribution is linear in dT (see also [12]), 

g W = 47TdT (34) 

whereas the remainder g^ cs ^ is bound by a constant independent of dT, both derived in detail 
in Appendix lAl 

Finally, in the Polyakov loop correlator (fT9|) . using Eqs. ([29]h (f3Tj) . ([33]) . and (fM|) . we obtain an 
exponential decay at large distance d = |r — r'| 

P(r)P(r')) = lexpf-^ + const.) (35) 



2 ^ V 2T 2 
or equivalently a linear growth of the free energy 

F QQ( d ) = ad + const - ( 36 ) 
and read off the string tension 

* = (37) 
2 T 

Given the dependence on p, d and T in Eqs. ([T9]) . ([2~T]) . ([25]) . and ([3T]) . the coefficient of a 
term linear in d can only be of that form (also for dimensional reasons). The achievement of 
this part of our work was to analytically prove this confining behavior and to determine the 
proportionality factor. 

As a side result we find that the holonomy dependence has dropped out completely in the infinite 
volume (technically because w enters together with J + , see Eq. (j!8p . this contribution, however, 
vanishes in the infinite- volume limit). This is consistent with the fact that the average Polyakov 
loop in our model actually vanishes for all holonomies in the infinite- volume limit, which is not 
difficult to show. 

In other words, the disorder generated by long-range fields of dyons dominates the effect of the 
holonomy on the average Polyakov loop. We remind the reader that this finding is based on the 
same density of all kinds of dyons for all holonomies. Hence our model is valid only at maximally 
non-trivial holonomy, i.e. in the low temperature phase, whereas in the high-temperature phase 
modifications are expected that may reintroduce a holonomy dependence. 



S 



3.3 Fixing the physical scale 

With the string tension of Eq. (j37[) at hand, we can set the scale of our model. All analytical and 
later numerical calculations provide, of course, relations between dimensionless quantities. As 
we have done already, we can measure all lengths in units of the inverse temperature j3 = 1/T. 
For the string tension this means 

* P « f = l {fp f. (38) 



T 2 2 T 3 2 VP~ 1/3 / 2 

The ratio on the left hand side is known from lattice simulations. We have introduced a "packing 
fraction" fp of the dyon gas, since p~ 1 ^ represents the mean distance and /3 can be interpreted 
as being proportional to the core-size of corresponding non-Abelian dyons. 

We resort to lattice results on the SU(2) string tension and its temperature dependence in |26j . 
We parameterize these results (cf. Fig. 3 in that reference) by 

-co - ,a -y v^ryi (39) 



a(T = 0) V T C J V V T t 

but additionally require a{T)/a{T = 0)|t=o — 1> which amounts to B = 1/A — 1. We find 
A = 1.39 to describe the lattice data reasonably well. 



Using another lattice result, T c j \fa{T = 0) 0.71 [27], allows to rewrite Eq. (|39|) according to 
a(T) _ *(T = 0)fT e yj, T\^f i + B /_T\^\ (4Q) 




T 2 T 2 \T \ T r \ V T, 



Together with (|38p this formula relates the density of dyons respectively their packing fraction 
to the temperature ratio T/T c . Finally, physical units can be introduced using a(T = 0) = 
(440 MeV) 2 (as we already did in [22]) corresponding to T c = 312 MeV. 

Assuming that our dyon gas model provides the correct phenomenological value of the string 
tension we can tell how the density and the packing fraction have to behave as functions of the 
temperature below T c . For both the limits T — > T c and T — > the density p tends to zero. Its 
maximal value p ma x ~ 0.25T c cr(T = 0) is reached at T ~ 0.65 T c . In physical units we have 
Pmax ~ 2 fm~ 3 . The packing fraction fp for our model diverges for T — > and tends to zero 
for T — > T c . The latter behavior can be interpreted such that the diluteness assumption applies 
best near the phase transition, but becomes more and more violated for low temperatures. This 
problem, however, is well-known to occur also for the instanton liquid model (see e.g. [1]). 



3.4 Polyakov loop correlator at arbitrary separation and finite-volume effects 

In this subsection we numerically evaluate the Polyakov loop correlator from Eq. (|19p and 
correspondingly the integrals I± from Eq. (|20|) at arbitrary quark-antiquark separation d and 
arbitrary volume V (both finite and infinite). This allows to investigate finite- volume effects. 



9 



As a by-product we will confirm the linear behavior for infinite volume and large separations, 
Eqs. flMEU). 

To perform the numerical integration efficiently, we split the integrals in two regions, S = S%„, 
a ball of radius R < R, and its complement S = S RT — 5: 

- XMK^^HMK^k^))- (41) 

The integral over S can be solved numerically with standard methods, e.g. ordinary Monte Carlo 
sampling, because both the region of integration and the integrand are finite. By introducing 
spherical coordinates it can even be reduced to a 2-dimensional integral: 

f±,S = [ ds expf U — — — -± 



s 

RT 

2 



2\|rT-s| |r'T-s| y 

1 1 



Jo 



2tt ds s I dQ sin 9 exp I — 



2\D + {s,9,dT) D^(s,6,dT) 



with D± according to Eq. (|22p . For R—too the integrals over S exhibit infinities, which need to 
be subtracted, before a numerical treatment is possible. For finite but large R this subtraction 
is essential for an efficient computation of the integrals. To exhibit the infinities, we expand in 
powers of 1/s: 

pRT p7T 

= 2tt dss 2 I dQ sin^ 
Jrt Jo 

1 + l_^ + < (-3W' + »WW,-4) +0(1/<4) 

V{S) + A + finite (42) 

ds C o S ( * _ _i -) = f ds cos (-( + 0(l/s 3 - 

§ V2|rT-s| 2\r'T-s\J J- s \2s\ s 2 w ' 

2tt J dss 2 J d9 sine ^l + C(l/s 4 )^ = y(5) + finite, (43) 



where 



RT r"K 



A = 2vr /_ dss 2 I de S me\ Z ~-^ + ^ ^l^;> — - (44) 



1 i (-3 {dT) 2 + 9 {dT) 2 cos 2 0- 4) 



AT J 

27T / „o / ~ o „n\ „ / ~ „\ / i? 



3^^-^ +3T^-i?J +iln - (45) 



see also (f26|h Note that the imaginary part of /_ vanishes, as argued in Section [3? 
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The finite parts of the above integrals can be evaluated numerically: 



-,5, finite 



-,5, finite 



[•RT f'TT 

2tt / dss 2 / dO s'm6 
Jrt Jo 



exp - 



+ 



_ i J_ 

s + 2s 2 

2vr / dss 2 dO sin(9 
Jrt Jo 



2\D + (s,e,dT) D^(s,9,dT) 
3 (dT) 2 + 9 (dT) 2 cos 2 9 - 4) 



(46) 



cos 



1 



2D + (s,6,dT) 2D_(s,9,dT) 



The range of integration of J ds, which extends to infinity in the limit R — > oo, still poses a 
problem, but can be overcome by a change of variables according to 



ds ■ 



dx 



(47) 



(we have chosen that particular form, because the integrands of /, § finite and /_ § finite are 
proportional to 1/s 2 for large s.) Consequently, 



/ (is' —p5=\ dx' 
Jrt s u J Xo 



1 / RT — x + xq 

For simplicity and without loss of generality we choose xq = in the following. Then 



(48) 



RT 



RT 



dsF(s) 



1/RT-l/RT F{ l/{l/R T - x )) 

dx = 

(1/RT-x) 2 



(49) 



where the integrand is roughly equally distributed over the finite range of integration < x < 
1/RT — 1/ RT, if F(s) ~ #/s 2 - The final expressions for numerical evaluation are 



,S, finite 



rl/RT-l/RT pir 

2ir / dxs A / d9 sin6l 

Jo Jo 



cxp 



1 



+ 



1 



i 1 
" 1_ s + 2? 



2\D + (s,9,dT) D^(s,6,dT) 
i (-3 (eff) 2 + 9 (dT) 2 cos 2 0-4) 



24s 3 



(50) 



-,5, finite 



rl/RT-l/RT rw 

2ir / dxs A d9 sin 6 

Jo Jo 



COS 



1 



2D + (s,9,dT) 2D_(s,e,dT) 



where s = 1/(1/ RT — x). 
In total 



/+ = U, s + f + ,sM^ + V-V(S)+A 
f- = f-f + f-awto + V-ViS). 



(51) 
(52) 
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The Polyakov loop correlator for maximally non-trivial holonomy is 



>(r)P(r')> = iexp^tflni^I) - \ exp Uk In J^H (53) 



(cf. Eq. (|17p ). In the limit V — > oo this equation simplifies to 

P{t)P(t)J = 2 exp l T3 )• ( 54 ) 

We have performed the remaining integrations numerically and show the results below in Figs. 
13 and O 



4 Ewald's summation method 
4.1 Outline of the method 

In the following we briefly summarize Ewald's method. For a more detailed presentation we 
refer to [28]. Our main motivation to use this method is to systematically control finite- volume 
effects in observables, in particular those contributing to the Polyakov loop in Eq. (|8j). 

The first step in Ewald's method is to mimic the infinite space by sampling the physical system 
restricted to a basic cell, the so-called "super cell" , of spatial volume L 3 which will - for finite 
density - contain only a finite number of randomly placed dyons. In a second step the space is 
filled with replicas of the super cell shifted by nL, n G Z 3 . Sums over infinitely many dyons in 
infinite space are replaced by sums over these replicas. 

The infinite sum <I> in the Polyakov loop, Eq. ([8]), is modified tod 

= EE, S — FT. ( 55 ) 

nGZ 3 j 1 J 1 

where j = (i, m) is now a superindex running over all dyons and antidyons coming in equal 
number (J takes no = 2K different values). 

Naively one might think that such a sum can be approximated by summing over a large but 
finite number of copies of the super cell. One can show, however, that even though this sum 
converges, when increasing the total volume further and further, it converges to a result that 
differs from the desired infinite sum ( I > (r) (cf. appendix [B]) . The distortion depends on details 
of the charge distribution like surface charges. Only in the limit L — > oo it is expected to be 
identical to the Ewald result. 

The third step and key idea of Ewald's method is to split the terms l/\r — tj — nL| in $ in a very 
specific way into an exponentially decaying "short-range part" and a smooth "long-range part" . 
While the sum over the terms appearing in the short-range part converges exponentially, the 



Note that Ewald's method is quite general in a sense that it is capable of performing infinite sums of arbitrary 
inverse powers |29] . 
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sum in the long-range part is carried out in Fourier space, where it also converges exponentially. 
This allows a rather efficient computation of the sum in Eq. (|55p up to arbitrary precision. 



In detail the splitting into the short and long-range sum is done in the following way: 

$(r) = $ short (r) +$ long (r) (56) 



y/2\ J J \r-Tj- nL\ 
-«(,) . EE e rf (t^)^^, (58) 



1GZ3 J 



where erf denotes the error function. The physical intuition behind this decomposition becomes 
clear by computing the charge corresponding to this potential, i.e. by applying the Laplace 
operator to 4>. Of course, the original potential l/|r — rj — rvL\ yields pointlike sources at the 
dyon positions rj + nL. The auxiliary term — erf(|r — — nL\/y/2X)/\r — r\- — nL\ yields a 
continuous charge distribution around the same locations, but with Gaussian profile of width 
A and opposite sign. It is clear that the effect of these two charge distributions increasingly 
cancels in <J> short with growing distance, actually in an exponential manner. 

In $ lon s the smeared charge generates a non-singular potential at the dyon positions. This leads 
to a convergence in its Fourier transform, which is exponential, too: 



$ lon g(r) = g £ e+ , k (n )r e " A k 2 ^ 2/2 (f>e-* k(n)rj ) , k(n) = ^n. (59) 



The expressions in parentheses are called structure functions, since they contain the information 
about the dyon positions. Note that this expression for $ long is correct only if the system is 
neutral, i.e. if ^ • qj = 0. This is the case for the non-interacting dyon model. For non-neutral 
systems $ long diverges. The free parameter A determines the tradeoff between the long-range 
sum and the short-range sum. While the short-range sum can be evaluated rather quickly for 
small A, the opposite is the case for the long-range sum. The optimal choice for A is discussed 
in the following section. 



4.2 Performance of Ewald's method 

To determine the free energy of a static quark-antiquark pair within the non-interacting dyon 
model, we need to evaluate Polyakov loop correlators. Doing this in an efficient way amounts 
to computing $ at a set of sample points r distributed on a cubic lattice throughout the spatial 
volume. Let M be the number of sample points. The computational costs to evaluate the 
short-range sum (|57p up to any desired accuracy are 0(MA 3 ) assuming that dyons within a 
spherical region around a given sample point can be always identified within the same cpu-time 
(see below, how this can be realized). 

Similarly one can read off the computational costs for evaluating the long-range sum ([59]) up 
to exponential precision. The structure functions ^ ■ qje~ lk ^ r ^> are independent of the sample 
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point r and, therefore, need to be computed only once for a given dyon configuration. The 
number of the required structure functions is proportional to V/X 3 , hence the corresponding 
computational costs are of order 0(V 2 / X 3 ). The time needed for the subsequent computation of 
3> long at all M sample points is proportional to MV/X 3 . Consequently, the total computational 
costs to perform the long-range sum are 0{V 2 /X 3 ) + 0(MV / X 3 ). 

The computational costs of the short-range sum and of the long-range sum depend on A in just 
the opposite way (as expected). One should choose A in an optimal way, such that the total 
computational costs are minimized. Obviously the optimal choice for A also depends on M. 
Since typically M oc V, as it is the case for our computations, the optimal value for A should 
be chosen according to A 3 oc y/V . Then the performance of Ewald's method is 0(V 3 ^ 2 ). This 
behavior has been confirmed numerically, cf. Fig. [TJ left panel. 



4 r 




1000 10000 100000 0.8 1.2 1.6 2 2.4 

V A 

Figure 1: Performance of Ewald's method. Left: a log-log plot of the computing time needed 
to evaluate the potential at M oc V sample points per dyon configuration as a function of the 
spatial volume V. The density of dyons is p = 1.0. The vertical axis is labeled such that one 
unit of cpu-time was needed to perform the computation for V = (10. 0) 3 corresponding to a 
number of dyons no = 1000. A was chosen according to A 3 oc VV. The straight line with slope 
3/2 illustrates that for large spatial volumes/dyon numbers Ewald's method indeed exhibits the 
expected 0(V 3 ^ 2 ) scaling. Right: the computing time as a function of the parameter A (in 
units of the inverse temperature) for i%d = 8000 and V = (20. 0) 3 . The vertical axis is labeled 
such, that one unit of cpu-time was needed to perform the computation at the optimal value 
A pt ~ 1-5. 

Of course, A 3 oc y/V is only a statement on how to increase A, when enlarging the spatial volume 
V. How to choose A for a given V such that the corresponding computing time is minimized, 
has to be determined by numerical experiment. In the right panel of Fig. Q] we show in an 
exemplary plot corresponding to no = 8000 and V = (20. 0) 3 the computing time needed to 
calculate the dyon potential $ as a function of A. Obviously, there is an optimal choice for A. 

Note that in the literature there also exists another version of the just explained "classical 
Ewald method", the so-called "particle mesh Ewald method" (cf. e.g. [29]). This version is more 
efficient, when the interaction energy of a system of positive and negative charges needs to be 
computed. However, for our problem at hand, the computation of the temporal gauge field 
there is no advantage with respect to performance. Since it is significantly simpler to implement, 
we resort to the classical Ewald method. 

For an efficient computation of the short-range sum $ short (r) it is mandatory to determine, 
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which dyons are located in a spherical region of given radius R centered around r in 0(1) 
computer time. To this end we divide the supercell into a grid of cubic subcells and generate 
for each subcell a list of the contained dyons. In addition we have implemented a function that 
determines all subcells, which are inside or which intersect the surface of the above mentioned 
ball. Then we call all those subcells for the dyons they contain. In this way we do not need to 
inspect all the dyons in the supercell and check whether their distance to r is smaller than R. 
Of course, the grid of subcells has to be sufficiently fine-grained, to be able to mimic a ball of 
radius R with rather small cubes. 

5 Numerical results 

5.1 Extracting the infinite volume string tension using Ewald's method 

We compute the free energy of a static quark-antiquark pair as a function of their separation 
from Polyakov loop correlators as described in Section [2j We keep the dyon density p and the 
temperature T fixed and perform computations for various dyon numbers no, corresponding to 
various spatial volumes V = rip/p of the super cell. The superposition of dyon potentials <& is 
calculated by means of the Ewald method as explained in Section HI We restrict ourselves to 
maximally non-trivial holonomy. 

Of course, the resulting free energies are different for different dyon numbers, i.e. spatial volumes 
of the super cell, because of finite- volume effects. In particular the dyon potential <]? is L-periodic 
along the three spatial directions, which obviously implies periodic Polyakov loops and loop 
correlators. Therefore, L has to be chosen sufficiently large to ensure that the free energy can 
be determined for quark-antiquark separations of phenomenological interest, typically a few fm, 
without being significantly distorted due to periodicity. 

In Fig. [2] we show quark-antiquark free energies for p/T 3 = 1.0 and dyon numbers 1000 < no < 
125000 (corresponding to 10.0 < Lp 1 / 3 < 50.0) as functions of the quark-antiquark separation 
dp 1 / 3 . We also show the analytic results for finite and infinite volume in this plot. Note that 
we express lengths in units of p 1 / 3 , which is the average dyon separation in a random dyon gas. 
The number of dyon configurations used for each dyon number ri£> is listed in Tabled! It can be 
seen that the free energies converge, when increasing ri£> (or equivalently Lp 1 / 3 ). This allows an 



n D 


LpV3 


configurations 


1000 


10 


1600 


8000 


20 


800 


27000 


30 


120 


64000 


40 


90 


125000 


50 


60 



Table 1: Number of random dyon configurations used for every simulation at fixed dyon number 
ri£) or dimensionless length of the volume in which the dyon positions are sampled, Lp 1 / 3 , 
respectively. 
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Figure 2: Free energy of a static quark- antiquark pair as a function of its separation for p/T 3 = 
1.0 and various supercell extensions Lp 1 / 3 corresponding do different dyon numbers ud- In 
addition we show the results obtained from a numerical evaluation of the analytic result at finite 
and infinite volume. For better visibility the analytic results are shifted by log 2 and therefore 
the corresponding curves start close to the origin. 



extrapolation to infinite volume. In the left panel of Fig. [3] we show linear extrapolations of the 
finite-volume static free energy to infinite volume for a number of quark-antiquark separations. 
We also compare the results of the extrapolation to the analytically obtained free energy at 
infinite volume in the right panel of Fig. [3l As can be seen, analytic and extrapolated results 
nicely agree within errors. 

Let us point out that there are other methods of obtaining an infinite volume result numerically 
without employing Ewald's summation method. 

An obvious method is a straightforward superposition of dyon potentials in a finite cubic box of 
size L 3 , that we call dyon sampling volume. We have used this method in a previous publication 
[22], to which we refer for further details. Note that there is no exact translational invariance 
anymore, in contrast to when using periodic boundary conditions via Ewald's method. To keep 
finite-volume effects at a tolerable level, we have to restrict the evaluation of Polyakov loops 
to a spatial subvolume sufficiently far away from the boundary of the dyon volume. We will 
call this subvolume field sampling volume. It is centered inside the dyon sampling volume and 
has extension I < L. On the one hand, finite- volume effects are expected to be negligible for 
sufficiently small I. On the other hand, however, decreasing I reduces the available information 
per dyon configuration and, therefore, reduces statistical accuracy. In practice one would need 
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Figure 3: Demonstration of the infinite volume limit. Left: Static free energy for ten distances 
as a function of the inverse supercell extensions Lp 1 / 3 (corresponding to different dyon numbers 
ri£)) and its extrapolation to infinite volume. Right: The latter compared to the analytic result, 
Eq. (|54p . for arbitrary distances in infinite volume. 



to identify plateaus in the observables as functions of I. An obvious and major drawback of 
proceeding in such a way is that one needs to extrapolate in two parameters, the extension / of 
the sampling volume and the extension L of the dyon volume. Clearly this is technically more 
complicated, than what has to be done using Ewald's method, where the only parameter subject 
to extrapolation is the extension of the supercell L. 

One can also think of evaluating just one Polyakov loop correlator in the center of the volume 
of each random configuration. We should point out that this is not really feasible if there are 
interactions, since a significantly larger statistics is needed when no volume averaging is done. 
For the non-interacting case it is applicable and therefore worth being mentioned. 



6 Summary and outlook 



In this work we have shown analytically that a non-interacting random dyon gas leads to a 
linearly rising free energy of a static quark-antiquark pair as a function of the distance in 
between. Correspondingly the string tension a turned out proportional to the ratio of the 
density and the temperature, i.e. to p/T, cf. Eq. ([37]) . We were able to present explicit formulae 
for arbitrary distance and for finite volume with certain integrals left to be evaluated numerically. 
We convinced ourselves that the dependence on the holonomy drops out in the infinite volume 
limit. This reflects the fact that - concerning the Polyakov loop and its correlator - the model is 
able to describe only the confinement phase. For the deconfinement transition as well as for the 
deconfinement phase, where the center symmetry becomes broken, the model should be altered 
taking into account that dyons with opposite charge should be statistically weighted differently. 

We emphasize, that our analytical approach is specific for the non-interacting case. For the 
interacting case it will not be applicable without approximations, and in the first instance we 
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will have to rely on numerical simulations. Strong finite-size effects of the naive treatment with 
finite boxes containing the dyon sources have led us to employ a numerical method well-known 
in the physics of a three-dimensional Coulomb plasma, the Ewald summation method. We 
convinced ourselves that this method will be applicable also to the more realistic interacting 
dyon gas. 

Indeed, we have demonstrated, how Ewald's summation method can be used to deal with long- 
range objects also in field theory, in our case with random ensembles of dyon constituents. 
In this semiclassically motivated model we have computed the local Polyakov loop, and from 
its correlator we have extracted the string tension, the main observable characterizing confine- 
ment /deconfinement at finite temperature. 

The Polyakov loop is a function of an infinite sum of Coulomb contributions of dyons with both 
signs of charge (cf. Eq. ([8])). According to Ewald's method we have decomposed this sum into 
short-range and long-range parts, Eqs. (157|) and (159p . and have optimized the width A of the 
auxiliary Gaussian charge cloud governing the strength of the exponential convergence of both 
parts. 

Figs. [2] and [3] show our main results, the free energy of a quark-antiquark pair as a function of 
its separation, for various extensions of the (periodically repeated) supercell volume, but fixed 
dyon density. These figures also demonstrate that the straightforward extrapolation to infinite 
supercell volume is a valid procedure to obtain results for an infinite non-interacting system of 
dyons. In this limit the Polyakov loop correlator behaves as expected: it decays exponentially 
toward larger quark-antiquark separations. The corresponding string tension can be read off 
unambiguously (and used to fix the physical scale of this model). 

We have discussed the advantages of Ewald's periodic summation over methods that at finite 
volumes measure observables only in subvolumes: it keeps translational invariance and the 
infinite volume limit amounts to extrapolating just one parameter. 

The applicability of the numerical method we have used is not restricted to a non-interacting 
dyon ensemble and/or to SU(2). Dyon fields in higher gauge groups decay with the distance in 
the same Coulombic manner, just possessing different color structures. Several other ingredients 
of dyon models contain Coulomb tails, too, like the interaction of dyons via the action or 
their moduli space metric. Furthermore, spatial Wilson loops (providing an area law decay with 
magnetic screening persistent also in the deconfined phase) can - with the help of Stokes' theorem 
and based on (anti)selfduality - be represented as area integrals over the normal component of 
the gradient of the same infinite sum. 

The ability to perform a controlled infinite volume extrapolation (with a single remaining pa- 
rameter L, the extension of the periodically continued spatial volume) is even more important 
in more complicated systems. An ensemble of random dyons could be treated easily with up 
to 10 5 dyons. Interacting dyon ensembles are numerically much more expensive such that the 
reduction to a smaller number of dyons most likely cannot be avoided. Then finite- volume effects 
might become a limiting factor. Consequently, Ewald's summation method seems to become 
indispensable, however, in form of the particle mesh Ewald method, which is more efficient than 
the classical Ewald method, when computing dyon interactions. 

Finally, one could think about applying Ewald's method to more complicated objects in gauge 
theory, whose corresponding fields are also of long-range nature, such as merons or regular gauge 
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instantons (30J E] and generalizations thereof [32j [33] . 



A Calculation of some integrals 

We derive the following results for the integrals and g^ cs ^ involved in Polyakov loop corre- 



lators in Section 13.21 

2 



[ ds ( - 1 , - , / \ = 4vr|u - u'l (60) 

7r3 VI u - s I l u - s l/ 

0< / ds \ cos ( , 1 , - — - — r ) - l + -( , , -— - — r ) 1< const. (61) 

7 R3 t Vl u - S l l u - S \J 2 VI U - S I l u — s|y J 

In the first integral we use the well-known Fourier representation of the Coulomb potential 

1 1 f 47T ■ 

dp^e tps (62) 



|s| (27I") 3 J p 
to calculate 

2 



1 1 



1 f , , 1 

4^4 



= ^y"dp^2-2cos(p(u-u')) J =8^ ^ d0sin0^1 - cos (p|u - u'| cos 0) J 

_ 8 rA( 2 - 2 ^frl"-"'l)) = 4 T |!-u? " (63) 
Jo P V p|u-u'| ) 

The Laplace operator can be used to check this result. Acting with respect to u' and u on the 
left hand side we obtain (from the mixed term in the integrand) 

A H A M , Jds { — 1 ^ - — - ^ = -2 Jds (-4vr) 2 5(u - s)5(u' - s) 

= -32vr 2 5(u - u') (64) 

On the right hand side it gives the same since 

A U A U / 4tt|u - u'| = A u 8?r = -32tt 2 (5(u - u) . (65) 



The integrand of the second integral cosx — 1 + x 2 /2 = h(x) is positive, which proves the first 
inequality. For the second inequality we split s-space into two half-spaces, |u — s| «s | u ' — s |> 
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separated by the midplane between the two points u and u'. The integral over each half-space 
gives half of the full integral and thus we can specify to one of them, e.g. where 



111 

< i r - r~, r < i r 66 ) 

|u — s| |u' — s\ |u — s| 

holds. Since the integrand h(x) is monotonically increasing for x > 0, we obtain an upper bound 

1 1 \ 1/1 1 

ds { cos , r - — r - 1 + - 



u — s| |u' — s| / 2\|u — s| |u' — s w 

1/1 1 



2 / ds { cos ( 1 r - — ; r ) - 1 + - 

u— s|<|u'— s| 



u — s| \u' — s\ J 2\|u — s| |u' — s 

< 2 / ds cos — *— - 1 + - 1 1 (67) 

u — s 2 u — sr 



Due to the positivity of the integrand, we can extend the latter integral back to the full space 
and by virtue of translational invariance put u = obtaining another bound 

ds <! cos I ! — - — r - ~i — t~ — r I - 1 + r ( 



f ds\ 



2 V u — s u' — s 



,f 1 11 2ir 2 

< 2 • 47r / ds s \ cos 1 + 







s Is 2 ) 3 

independently of |u — u'|, by which we have proven the second inequality. 



B Ewald's sum compared with summing over a finite array of 
supercells 

Ewald's method amounts to summing over infinitely many copies of the cubic spatial volume L 3 
called supercell. An alternative approach is to truncate this sum after a large but finite number 
of copies in every spatial direction dbx, =by and ±z. The corresponding dyon potential obtained 
by summing over (2re + l) 3 copies of the supercell is then 

+ra +n +n 

* flnitesum w = E E E Eu-£ 



r, — nL 

n x =—nn y =—nn z =—n j ' J ' 

where r, Tj € [— L/2, +L/2] 3 and n = (n x ,n y ,n z ). One might expect that, when n is chosen 
sufficiently large, the Ewald result, denoted by <J> Ewald i and $ fimtc sum become arbitrarily close. 
In this section we explain that this is not the case, i.e. that even though $ fimte sum converges, 
when increasing n, it will in general differ from $ Ewald . 

The difference between the two approaches is 
A$(r) = $ Ewald (r) - $ finite sum ( r ) 



20 



E 

j{— n,...,+n} n y &L 



E 

n,...,+n} n z G2 



E E 

,{-n,...,+n} j 



Vj - nL\ 



(70) 



In the following we demonstrate by means of a simple example that A<J>(r) 7^ in general. To 
this end consider no = 2, a dyon (q± = +1) at position ri = (— ri/2,0,0) and an antidyon 
(f/2 = —1) at position Y2 = (+d/2,0,0). 

For d = L dyons and antidyons in (j70|) cancel exactly with exception of antidyons/dyons located 
on planes at (^(n + 1/2), n y ,n z )L, n y ,n z £ Z. Since the dyon potential is identical to the 
potential of an electric charge in classical electrostatics, the situation is reminiscent to that 
of a uniformly polarized cubic dielectric with volume {{2n + l)L) 3 . For n ^> 1 the discrete 
charges can be approximated by the surface charge density a = ±4-7r/L 2 at the two opposite 
sides x = ±(n + 1/2)L. 

For d < L the dyon and antidyon potentials only partly cancel resulting in a reduced surface 
charge density a = ±47rd/L 3 . 

For n>l the difference A$ is given by 



r+(n+l/2)L 

V(A$(r)) = -2 / dy 

J-(n+l/2)L 



+ (n+l/2)L 
(n+l/2)L 



dz 



d(r-((n+l/2)L,y,z) 
4vrL 3 |r- ((n + l/2)L,y,z) 



(71) 



i.e. 



A$(r) 



4ird 
3L3 : 



x(l + 0(l/n] 



(72) 



In Fig. H] we show that this analytical result is accurately reproduced by our numerical imple- 
mentation of Ewald's method and the finite sum (|69|) using n = 50. 
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Figure 4: LA<I> as a function of x/L [y = z = 0) for a dyon at (—d/2, 0, 0) and an antidyon at 
(+d/2,0,0) (cf. text for details), a) d = L. b) d = L/2. 



For a larger number of dyons with arbitrary positions A$ is, of course, rather hard to estimate 
analytically. The physical picture, however, will remain the same: like in a polarized dielectric 
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surface charges will cause a difference between <j) Jwald an d $ hmtc sum . Only in the limit nn — > oo 
corresponding to L — > oo both approaches are expected to become identical. 

In principle both approaches can be used to simulate dyon ensembles, since, after appropriately 
extrapolating the dyon number nu — > oo (or alternatively L — > oo), one should obtain the same 
correct infinite volume result. We consider, however, Ewald's method to be superior, because in 
this approach the spatial volume is translationally invariant. This allows to maximally exploit 
a given dyon gauge field configuration by evaluating observables throughout the whole spatial 
volume. In contrast to that, translational invariance is broken when truncating the sum over 
copies of the super cell. Observables must only be evaluated in regions, where this breaking is 
sufficiently mild. Each observable requires to determine a corresponding region of sufficiently 
mild finite volume effects. Moreover, one has to assure that the associated systematic is removed 
by the infinite volume extrapolation. 
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